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ABSTRACT 

We present new .ffersc/id-SPIRE imaging spectroscopy (194-671 fim) of the bright starburst galaxy 
M82. Covering the CO ladder from J = 4 — > 3 to J = 13 — > 12, spectra were obtained at multiple 
positions for a fully sampled ~ 3 x 3 arcminute map, including a longer exposure at the central 
position. We present measurements of 12 CO, 13 CO, [C I], [N n], HCN, and HCO + in emission, along 
with OH + , H20 + and HF in absorption and H2O in both emission and absorption, with discussion. 
We use a radiative transfer code and Bayesian likelihood analysis to model the temperature, density, 
column density, and filling factor of multiple components of molecular gas traced by 12 CO and 13 CO, 
adding further evidence to the high- J lines tracing a much warmer (~ 500 K), less massive component 
than the low-J lines. The addition of 13 CO (and [C I]) is new and indicates that [C I] may be tracing 
different gas than 12 CO. No temperature/density gradients can be inferred from the map, indicating 
that the single-pointing spectrum is descriptive of the bulk properties of the galaxy. At such a high 
temperature, cooling is dominated by molecular hydrogen. Photon-dominated region (PDR) models 
require higher densities than those indicated by our Bayesian likelihood analysis in order to explain 
the high-J CO line ratios, though cosmic-ray enhanced PDR models can do a better job reproducing 
the emission at lower densities. Shocks and turbulent heating are likely required to explain the bright 
high-J emission. 



1. INTRODUCTION 

M82 is a nearly edge-on galaxy, notable for its spectac- 
ular bi£c4a£_out^wjjid high IR luminosity (5.6 x 10 10 
L , iSanders et al.ll2003D . Though its high inclination 
angle of 77° makes it difficult to determine, M82 is 
likely a SBc barred spiral galaxy with two trailing arms 
(jMayya et al.l 12005) . Its r edshift-independent di stance 
is about 3.4 ± 0.2 Mpc (jDalcanton et al.l 12009ft . and 
after correcting the commo nly cited redshift (0.000677, 
Ide Vaucouleurs et al.l (|1991h ) with WMAP-7 parameters 
to the 3K CMB reference frame, we find a redshift of 
0.000939. Given this distance we assume a conversion of 
17 pc/". 

Due to its proximity, M82 is an exceptionally well- 
studied starburst galaxy. High star formation rates (9.8 
Mpt yr" 1 , likely enhanced by interactions with M81, 
lYun et al.j[l99l and a large gas reservoir produce bright 
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molecular and atomic emission lines. Such lines can yield 
important information on the interaction between the 
interstellar medium (ISM) and star formation (SF) pro- 
cesses, such as the influence of SF on the ISM through 
photon-dominated region (PDR) or other excitation, as 
traced by intermediate to high-J CO rotational lines. 

Ground-based s t udies of CO in M82 are numerous 
(IWild et all [1991: IMao et al.l [20001 IWeifi et al.1 [200l 
iWard et all 120031 : IWeifi et al l 120051 ) . examining both 
morphology and physical conditions of the gas. High- 
resolution CO maps of the 1 kpc disk indicate that the 
emission is largely concentrated in three areas: a north- 
east lobe, southwest lo be, and to a lesser extent, a central 
region (see Figure 1 of IWeifi et al1l200"l . The two lobes 
are separated dynamically, a s can be seen in position- 
velocity diagrams (Figure 3 of lWeifi et al.l[200l . Outside 
of the disk, m olecular gas emission is also detected in the 
halo/outflow (jTavlor et alj|2001l : IWalter et al.ll2002D . 

In addition to examining the morphology of molecu- 
lar gas, CO emission lines can be used to determine the 
physical conditions of the molecular gas in galaxies. Pre- 
viously, due to terrestrial atmospheric opacity, only the 
first few lines in the CO emission ladder could be studied. 
The first studies of higher- J lines (described below) have 
indicated that they can trace components of gas sepa- 
rate from those measurable with low-J lines. Many of 
the most interesting questions about galaxy formation, 
evolution, and star formation concern the balance of dif- 
ferent energy sources, i.e. what role might cosmic rays, 
ultraviolet light from stars, X-rays from powerful AGN, 
or turbulent motion play in the star formation history of 
various galaxies? In what way does star formation influ- 
ence the molecular gas, and vice versa? Estimating the 
influence of these various energy sources, however, often 
depends on knowing the physical conditions of the gas. 



2 



KAMENETZKY ET AL. 



We therefore model physical conditions of these high-J 
lines first in order to inform our later discussion on energy 
sources. Other molecules are also useful in this study; in 
a ground-based sur vey of 18 different molecular species, 
lAladro et all (|2011[ ) also found that some molecules trace 
different temperature components than others and that 
the different chemical abundances in M82 and NGC253 
may indicate different evolutionary stages of starbursts. 

The Herschel Space Observatory (|Pilbratt et al.|[20Tol) 
is the unique facility that can measure the submillime- 
ter properties of nearby galaxies in a frequency range 
that cannot be observed from the ground. As one of the 
brightest extragalactic submillimeter sources, M82 has 
been studied extensively with Herschel. For example, 
the imaging photometer of the Spectral and Ph otometric 
Imaging REceiver (SPIRE. iGriffin et al.ll20TfH ) has been 
used to study the cool dust of M82, revealing wind/halo 
temperatures that decrease with distance from the cen- 
ter wi th warmer starburs t-like filaments between dust 
spurs (Roussc l et al.ll2010[) . The tidal interaction with 
M81 was likely very effective in removing cold interstellar 
dust from the disk; more than two thirds of the ex t rapla - 
nar dust follows the tidal streams. iPanuzzo et al.l (|2010f ) 
used the SPIRE Fourier- Transform Spectrometer (FTS) 
to study a single spectrum of the 12 CO emission from 
J = 4 -> 3 to J = 13 -> 12 to find that these higher-J CO 
lines likely trace a ~ 500 K gas component not seen in the 
~ 30 K component that can be observed from ground- 
based studies. Also on-board Herschel is the Heterodyne 
Instru ment for the Far Infrared (HIFI, Ide Graauw et al.l 
(2010)), which consists of a set of 7 heterodyne receivers 
with resolution of 125 kHz to 1 MHz for electronically 
tuneable frequency coverage of 2 x 4 GHz; it covers 480 
- 1910 GHz. HIFI found ioni zed water absorption from 
diffuse gas (|Weifi et al.ll2010D and high-J transitions of 
the CO ladder. These CO transitions indicated a combi- 
nation of one low and two high density gas compon ents 
via comparison to PDR models (jLoenen et al.ll2010f l. 

We confirm the presence of multiple molecular hydro- 
gen thermal components in M82 by performing a more 
in-depth modeling analysis on a deeper dataset as part 
of the Herschel Very Nearby Galaxies Survey. We add to 
existing data by using the SPIRE FTS mapping mode, 
providing spectroscopic imaging of a region approxi- 
mately 3' x 3', which helps us confirm our source-beam 
coupling corrections. We also pr esent a deep pointed 
spectrum (64 scans vs. 10 scans in IPanuzzo et al.l I2010T) 
in order to detect fainter lines, such as 13 CO, H2O, OH + , 
HF, and more. 

We add depth to the analysis by modeling both the 
cool and warm components of molecular gas, and si- 
multaneously accounting for 12 CO, 13 CO, and [C I]. 
We also use [C I] emission as a separate estimate 
of total hydrogen mass and other absorption lines for 
column density estimates. We first analyze the CO 
excitation using likelihood analysis to determine the 
physical conditions, and then compare to possible en- 
ergy sources. Our likelihoods test the uniqueness 
and un certainty in t he co n ditions, as has also bee n 
done inlNaylor et all (120101): iKamenetzky etHI (120111): 
iScott et all d201 ID; Bradford et al. l (|2009H ; |Panuzzo et al.l 
(120101k iRangwala et al.l (|201lD . 

Our observations are described in Section [2] We de- 



scribe the Bayesian likelihood analysis used to find the 
best physical properties of the molecular gas in Section [3] 
and present the results in Sections 14. II and !4. 21 In the re- 
mainder of Section [4j we discuss absorption results that 
are new to this study, possible excitation mechanisms of 
the warm gas, and comparisons to other galaxies. Con- 
clusions are presented in Section [5j 

2. OBSERVATIONS WITH SPIRE 
2.1. The SPIRE Spectrometer 

The SPIRE instrument (iGriffin et al.ll2010D is on-bo ard 
the Herschel Space Observatory (jPilbratt et al.ll2010l ) . It 
consists of a three-band imaging photometer (at 250, 350, 
and 500 /im) and an imaging Fourier-transform spec- 
trometer (FTS). We are presenting observations from the 
FTS, which operates in the range of 194-671 /im (447- 
1550 GHz). The bandwidth is split into two arrays of 
detectors: the Spectrometer Long Wave (SLW, 303-671 
/im) and the Spectrometer Short Wave (SSW, 194-313 
/im). The SPIRE spectrometer array consists of 7 (17) 
operational unvignetted bolometers for the SLW (SSW) 
detector, arranged in a hexagonal pattern. In the SLW, 
the beam FWHM is about 43" at its largest, dropping 
to 30" and then rising again to 35" at higher frequency. 
The SSW beam is consistently around 19". 

Two SPIRE FTS observations from Operational Day 
543 were utilized in this study: one long integration 
single pointing of 64 scans total ("deep", Observation 
ID 1342208389, 84 min [71 min integration time], AOR 
"SSpec-m82 -deep") and one fully-sampled map ("map", 
Observation ID 1342208388, - 5 hrs, AOR "Sspec- 
m82"). The map observation was conducted in high- 
resolution (HR) mode and the deep observation was con- 
ducted in high+low-resolution (H+LR) mode. Both were 
processed in high-resolution mode (Ai/ ~ 1.19 GHz) with 
the Herschel Interactive Processing Environment (HIPE) 
7.2.0 an d the version 7.0 SPIRE calibration deriv ed from 
Uranus (jSwinvard et alJl2010b iFulton et al.ll201fl ). 

2.2. Spectral Map Making Procedure 

In mapping mode, the SPIRE detector arrays are 
moved around the sky to 16 different jiggle positions, 
creating 112 and 272 spectra of 16 scans each for SLW 
and SSW, respectively, covering an area of the sky ap- 
proximately 3x3 arcminutes. The positions of these 
scans on the sky are presented in Figure [TJ with blue 
asterisks for SLW and red diamonds for SSW. 

The recommended map making method bins the spec- 
tra into pixels approximately one half the FWHM of the 
beam for each detector, which are about 35" for SLW and 
19" for SSW, leading to pixel sizes of 17.5"and 9.5". We 
wrote a custom script to create the map, based largely on 
the NaiveProjection method described in the SPIRE 
documentation for HIPE. Each of the 256 scans per de- 
tector were processed individually. All scans for a given 
detector and jiggle position falling within a given pixel 
were then averaged and an error bar for each wavenumbcr 
bin average is determined as the standard deviation of 
the scans divided by the square root of number of scans. 
All detector averaged spectra that fall into a pixel are 

11 SPIRE Data Reduction Guide, 

http: / /herschel. esac.esa.int /hcss-doc-7.0/ print / spire_dum/spire_dum.pdj 
December 2, 2U11. 
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Fig. 1. — Spectrometer Mapping Point Locations. SLW spectra 
locations are indicated by blue asterisks, SSW by red diamonds. 
The pixel boundaries, spaced 9.5" apart, are indicated by solid 
black lines. Ha contours are also plotted in black to indicate the 
orientation of the (nearly e dge-on) galactic dis k, from the Mount 
Laguna 40 inch telescope (Cheng ct al. 1997). Contours are in 
decreasing intervals of 0.2 log(maximum), i.e. 10°, 10 — , 10 — 4 
... 10~ ' 2 . The green circle indicates the size of the 12 CO J = 4— >3 
43" beam FWHM. The black X marks the position of the single 
pointing ("deep") spectrum. 

then averaged using a weighted mean (where the weight 
is the inverse of the square of the error bar). We used 
the same 9.5" grid for both bands. 

Using a 9.5" grid introduced more blank pixels in the 
SLW map, because the SLW map is more sparsely sam- 
pled because of its larger beam, as can be seen in Figure 
[TJ However, this enabled the comparison of the same 
spectral locations across both bands (where the data 
were available), without averaging together spatially dis- 
crepant spectra in the SSW, as would happen if pixels 
were made larger. We emphasize that the blank pix- 
els are somewhat artificial; the whole galaxy has been 
mapped, and the pixels locations are simply meant to 
indicate the central location of the detectors, though the 
beam size is larger than the pixel boundaries. 

2.3. Line Fitting and Convolution Procedure 

The mirror scan length determines the spectral resolu- 
tion of the spectrum. Because the scan length is neces- 
sarily finite, the Fourier- Transformed spectrum contains 
ringing; therefore, the instrument's line profile is a sine 
function, as can be clearly seen in Figures [3] and [31 The 
spectrum also contains the underlying continuum which 
must be removed before fitting the lines, which we do 
sequentially rather than simultaneously (see exceptions 
below). We isolate ± 25 GHz around the expected line 
center, and mask out the ± 6 GHz around the line cen- 
ter. We then fit the remaining signal with a second order 
polynomial fit to determine the continuum shape. After 
subtracting this continuum fit, we then use a Levenberg- 



Marquardt least-squares method to fit each line as a sine 
function with the following free parameters: central fre- 
quency, line width, maximum amplitude, and residual 
(flat) baseline value. The baseline value stays around 
zero because the continuum has already been subtracted. 
The central frequency is limited such that the line cen- 
ter is no greater than ± 300 kms -1 from the expected 
frequency given the nominal redshift of M82. For com- 
parison, the resolution varies from 230 to 810 kms -1 , 
from the shorter to longer wavelength ends of the band. 

For the deep spectrum, we detect weaker lines than 
in the map spectra. However, ringing from the strongest 
lines can interfere with the signal; therefore we first fit the 
strong lines ( 12 CO, [C I], and [N ii]) using the procedure 
outlined above and subtract their fitted line profiles from 
the spectrum. After all of the 12 CO, [C I], and [N n] lines 
are fitted and subtracted from the spectrum, we then do a 
second pass to fit the weaker lines. An illustration of the 
difference this process can make for the 13 CO J = 5^4 
line is in Figure [4] 

In general, all of the lines are fit independently, with a 
few exceptions: the 12 CO J = 7 -> 6 and [C I] J = 2 -> 1 
lines are a mere 2.7 GHz away in rest frequency, and 
ground state P-H2O and o-H20 + lines are separated 
by only 2 GHz. These two pairs must be fit simulta- 
neously. Both are fit independently to supply initial 
guesses, which are then used to fit both lines as the 
sum of two sine functions, each with their own central 
wavelength, width, and amplitude, but with one base- 
line value. 

The integrated flux is simply the area under the fitted 
sine function, which is proportional to the product of 
the amplitude and line width (converted to km/s). The 
error in the integrated flux is based on propagating the 
errors from the fitted parameters themselves. We note 
that the error estimation assumes all wavenumber bins 
are independent of one another, but that is in fact not 
entirely true in a FT spectrometer. Though lines that 
are separated spectrally do not affect one another greatly 
(hence why we fit most lines independently) , within each 
line fit the data points used in the 50 GHz range around 
the line center are not independent. 

The beam size of the SPIRE spectrometer varies be- 
tween the two detector arrays. In addition, it varies 
across the spectral range of the SLW, as described in 
Section f2.lt an d is not strictly proportional to wavenum- 
ber due t o the presence of multiple modes in the SLW 
detectors (|Chattopadhyay et alJ[2003T ) . When examining 
the spectral line energy distribution (SLED), it is imper- 
ative to scale all fluxes to a single beam size. For the map 
observation, we first fit the spectra as they were (with no 
correction factor). An example of integrated flux map, 
prior to any convolution or beam correction, is presented 
in Figure [21 with all other integrated flux maps available 
in the online version of the Journal. For the SLW detec- 
tor, we present integrated flux maps using both 9.5" and 
17.5" pixels. 

The integrated flux maps for each line were then con- 
volved with a kernel that matched the beams to the 
beam of the 12 CO J = 4 -> 3 map, which has a FWHM 
of 43". The kernel was created u sing a modifi e d ver - 
sion of t he procedure decrib ed by iBendo et al.l (|2011[ ) 
(see also iGordon et al.l [20081 ) . However, inst e ad of di- 
rectly applying Equation 3 from IBendo et al.l (|2011l ) to 
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Fig. 2. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 4— >3. Figures 2.1-2.19 are available in the online version 
of the Journal. The top half includes no beam correction or convolution. Black corresponds to the lowest flux or zero if any fluxes are 
negative, at which point the colorbar becomes purple. The bottom half is a map of signal/noise, though the color bar tops out at 20 in 
order to better illustrate which pixels are near the threshold of detectability. On the color bar, black corresponds to the lowest signal/noise 
or three if any pixels have S/N less than three, at which point the colorbar becomes purple. 
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the images of the beams, we applied the equation to one- 
dimensional slices of the beams to create the radial pro- 
file of convolution kernels, which we then used to create 
two-dimensional kernels. The approach worked very ef- 
fectively for creating kernels to match beams observed 
with SSW to the 12 CO J = 4— ^3 map. However, in cases 
where we created convolution kernels for matching beams 
measured at two different wavelengths by the SLW array, 
we needed to manually edit individual values in the ker- 
nel radial profiles to produce effective two dimensional 
kernels. 

After the mathematical convolution of the integrated 
flux map with the kernel (resampled to match our map 
sizes), the fluxes were all converted to the units of Jy 
km s _1 /beam, referring to the 12 CO J — 4 —> 3 beam. 
The ratio of beam areas was determined empirically by 
convolving the kernel with the smaller (observed) beam 
peak normalized to 1 and determining the ratio required 
to produce the larger Ocoj=4-^3 beam with the same peak 
(simulating the observation of a point source of 1 Jy km 
s^ 1 ). To account for blank pixels, only the portion of 
^kernel that encompasses data was used in the aforemen- 
tioned conversion. 

For the deep spectrum, w e used the same sourc e-beam 
coupling factor {r] c {v)) as in lPanuzzo et al.l (|2010f) . which 
was derived by convolving the M 82 SPIRE photometer 
250 /im map (|Roussel et al.|[2010h with appropriate pro- 
files to produce the continuum light distribution seen 
with the FTS. The deep spectrum was multiplied by 
this factor before the line fitting procedure. The deep 
spectrum is presented in four parts in Figure [31 and the 
measured lines fluxes (> 3c) are in Table Q] 

The central pixel of the convolved maps offers a direct 
comparison to the source-beam coupling corrected deep 
spectrum. In the SLW, these two SLEDS are within ± 
16% of one another. Later, we assume calibration error of 
20%, so these differences are within those bounds. There 
is greater variation in the SSW band, though this is the 
region in which the signal to noise of the lines greatly 
drops. When 20% calibration error is included, all the 
line measurements have overlapping error bars with the 
exception of 12 CO J = 13— >T2, where the deep spectrum 
measurement is more than twice that of the map. The 
12 CO J = 13 -> 12 map, however, has a S/N of only 4 
in the central pixel, with only 60/399 pixels having S/N 
greater than 3. We primarily use the deep spectrum for 
our likelihood analysis because of the higher S/N and 
access to more faint lines, but compare with using the 
convolved map central pixel in Section 14.21 The simi- 
larity of the two SLEDs, within error bars, using two 
independent methods (the derived rj c (u) from photom- 
etry comparisons vs. map convolution) to account for 
source-beam coupling, indicates that both methods are 
robust. 

Though the maps do not provide adequately high spa- 
tial/spectral resolution for a detailed study of the mor- 
phology of the emission, some qualitative assessments 
can be made. For example, the line centroids do trace 
the relative redshift /blueshift of the northeast and south- 
west compo nents {vhei ~ 300 k ms' 1 and 160 kms -1 , 
respectively. ILoenen et al.l [20101. However, we do not 
resolve the two separate peaks in flux. The capabilities 
of these maps to resolve gradients in the physical param- 



TABLE 1 

Measured Fluxes of Detected Lines in Deep 
Spectrum 



Transition 


^rest 

[GHz] 


[10 3 Jy kms" 1 ] 1 


12 CO 1 — 4— >3 


461.041 


85 66 -4- 90 
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O t U.iUO 
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12pn T — (\ v K 
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691.473 


7Q z.a -U n AA 


12pn T — 7 v « 


OUU.UJi 
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12 r , n 7 — « y 7 


Q91 son 
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1 n^fi Q1 9 


ao on -1- ft 7A 


J — 1U — t y 


llol.yoo 


on fin 1 n -1 1 
zy.ou ± U.41 


\^>\J J — 11 — > iu 


IzD i .U14 


iy.00 ± U.00 


12nr\ t in ,ii 

L;U J — LZ — ?-ll 


1 001 nnc 

looi.yyo 


-1 of I ri rr 
lo.zo ± U.oO 


12 f^r\ T 1 Q vIO 


14yo.yzo 


1 n a J_ 1 1 
1U.4 ± 1.1 


[CI] °Fi ->°ro 


492.161 


no no _L fi 

zo.SJo ± U.00 


[CI] i3 P 2 ^' 3 Pi 


809.342 


38.88 ± 0.58 


[Nil] 3 Pi-> 3 Po 


1462.000 


82.4 ± 1.2 


13 CO J = 5->4 


550.926 


3.83 ± 0.58 


13 CO J = 6^5 


661.067 


3.02 ± 0.15 


13 CO J = 7^6 


771.184 


1.84 ± 0.31 


13 CO J = 8->7 C 


881.273 


1.16 ± 0.45 


HCN J = 6^5 C 


531.716 


1.15 ± 0.42 


HCO+ J = 7^6 


624.208 


1.08 ± 0.17 


OH+ N = l->0, J = 0-»1 


909.159 


-5.2 ± 1.1 


OH+ N = l->0, J = 2->1 


971.805 


-8.88 ± 0.41 


OH+ N = l->0, J=l->1 


1033.118 


-9.94 ± 0.39 


HF J = 1^0 


1232.476 


-3.64 ± 0.34 


O-H2O 3l2 ->3o3 


1097.365 


1.39 ± 0.30 


0-H2O 3i2^2 2 i 


1153.127 


2.57 ± 0.37 


p-H 2 2iH-2 02 


752.033 


1.86 ± 0.36 


p-H 2 2 2->lll 


987.927 


3.09 ± 0.56 


P-H2O In ->Ooo 


1113.343 


-2.69 ± 0.44 


P-H2O 2 20 ^2ii 


1228.789 


2.01 ± 0.34 


0-H2O+ ln^Ooo 


1115.204 


-2.82 ± 0.42 



a All fits have been referenced to t he 4 3" beam of 12 CO 
J = 4 — > 3, as described in Section 12.31 Uncertainties do 
not include calibration error. 

b To convert to other units, we use these equations. L [Lq] 
= / F v dv [Jy kms" 1 ] X 0.012 u GHz . J Tdv [K kms" 1 ] 
= / F u dv [Jy kms" 1 ] 660.8 Vq\ z . F [W/m 2 ] = / F v dv 
[Jy kms" 1 ] X 3.3 xlO- 23 ^ GHz - 

c The 13 CO J = 8 -> 7 and HCN J = 6 -> 5 lines are detected 
at slightly less than S/N of 3; 2-<r upper limits would be 
9.0 and 8.4 x 10 2 Jy kms -1 , respectively. 

eters modeled in this work will be discussed in Section 
PI 

3. BAYESIAN LIKELIHOOD ANALYSIS 

We follow the method described in iWard et "all ([2003D 
for the 12 CO J = 6 -> 5 map of M82, used fre- 
quently in the analysis of ground based molecular ob- 
servations of galaxies (e.g. the Z-sp ec co llaboration 
i Navlor et all 120101: iKamenetzkv et al.ll2011t iScott et al.l 
120111 : [Bradfor d et al.1120091 ) an d also of the single p oint- 
ing SPI RE spectrum of M82 (|Panuzzo et all 120101 ) and 
Arp220 (|Rangwala et al.ll2011h . The goal of our Bayesian 
likelihood analysis is to map the relative probabilities of 
physical conditions over a large parameter space; this 
provides a more complete statistical analysis of the phys- 
ical conditions as opposed to simply finding one best-fit 
solution. 

For each molecular species, we first calculate a grid of 
expected line emission temperatures for various combina- 
tions of temperature (Tki n ), density (n(iJ 2 )), and column 
density per unit linewid th (Ni2 CO /dv) using RADEX 
(jvan der Tak et al.ll2007l ). We use the uniform sphere ap- 
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Fig. 3. — Deep Spectrum, split up by frequency. The switch from SLW to SSW bands occurs at approximately 950 GHz. Error bars are 
not shown on the spectrum for clarity, though the median error bar of each panel (times a factor of 20) is s how n in the upper left corner. 
The spectrum contains ringing because the line profile of the FTS is a sine function, as discussed in Section I2.3I Emission/absorption line 
locations are color-coded by molecular species: blue for 12 CO, red for 13 CO, light blue for water and its ion, orange for OH + , and violet 
for all others. 
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Fig. 4. — Example of line fit, using the 13 CO J = 5 — > 4 line. 
The top panel shows the continuum-subtracted spectrum (black); 
the ringing of the nearby 12 CO J = 5—^4 line is clearly visible on 
the right and interfering with the fit, overplotted in solid red. The 
middle panel shows the spectrum after subtracting out the strong 
12 CO, [C I], and [N n] lines (though only the subtraction of the 
12 CO J = 5—)- 4 line is visible in this figure). In this example, the 
signal to noise ratio of the fit was improved from 1.1 to 6.6 (where 
the signal is the integrated flux and the noise is based on propa- 
gating the statistical errors of the fit parameters in the integrated 
flux calculation). The bottom panel is the residual of the fit. In all 
panels, a horizontal line indicates zero. In the top two panels, the 
solid vertical line indicates the expected line center, and the two 
dashed vertical lines demarcate the area within that was not used 
for fitting the continuum. 

proximation for calculating the escape probabilities; the 
actual morphology in M82 shows a more complex velocity 
structure, therefore this approximation is considered an 
average of the bulk properties of the gas (and the results 
are not sensitive to uniform sphere vs. LVG approxima- 
tion). RADEX performs statistical equilibrium calcula- 
tions of the level populations, including the effects of ra- 
diative trapping, for a specified gas temperature, density, 
and column density per unit linewidth. The resulting so- 
lutions are output in the form of background-subtracted 
Rayleigh- Jeans equivalent line radiation temperatures. 

We use a 2.73 K blackbody to represent the cosmic mi- 
crowave background (CMB). We also experimented with 
using the continuum flux measured in our deep spectrum 
as a background; we fit the continuum in Jy across both 
bands, masking out the lines, with a third-order polyno- 
mial. The choice of this fit was to accurately represent 
the continuum background as a function of frequency; it 
was not meant to represent any physical conditions. Be- 
cause the relevant radiation background is the specific in- 
tensity (Jy/sr), we divide our continuum by the area cor- 
responding to a 43" beam (as the entire deep spectrum 
has been corrected to that) . Such a background is in fact 
orders of magnitude higher than the CMB at the highest 
frequencies that we are modeling. However, a grid with 
this background vs. the CMB produces the same likeli- 



TABLE 2 

RADEX Model Parameters and Ranges 



Parameter 



Range 



# of Points 



T kin [K] 


10 0.7 _ 1Q 3.7 


61 


n(i?2) [cm 3 ] 


10 2 - 10 6 


11 


Ni2 CO /AV [cm- 2 ] 


10 14 - 10 18 


11 




Kr 3 - io° 


41 


X ^co/ Xl2 CO 


ur 2 - lO" 1 


11 


X [CI]/ X ^CO 


lO" 2 - 10 2 


11 


AV [km s" 1 ] 


1.0 


fixed 


^background 


2.73 


fixed 



Note. — All parameters are sampled evenly in 
log space. Velocity is fixed because all modeling is 
column per unit linewidth. 

hood results. This is because at a kinetic temperature 
of ~ 500 K (the warm component we will model), colli- 
sional excitation greatly dominates over radiative excita- 
tion. In other words, at high temperatures, the modeled 
line intensities do not depend on the background radia- 
tion field. The cool component we will model is traced 
by low-J lines, whose background is not as affected, and 
so this component also finds the same results with cither 
background. Therefore we are presenting results using 
the CMB background. 

In addition to 12 CO, we also model 13 CO and [C I] as a 
function of the same temperatures, densities, and column 
density of 12 CO. For those molecular/atomic species, we 
also add the parameter of the relative abundance, e.g. 
[ 13 CO]/[ 12 CO] or Ai3 CO /Ai2 C0 . When modeling the 
intensity, the column density of 13 CO is simply that of 
12 CO times the relative abundance. [C I] is modeled 
with the parameter of [C I]/[ 12 CO]. Finally, we introduce 
one more parameter, the area fractional filling factor (&a- 
The modeled emission may not entirely fill the beam, 
so the flux may be reduced by this factor. All model 
grid points are therefore multiplied by each value of &a- 
The ranges of parameters, as well as the number of grid 
points, are presented in Table [2l 

The RADEX grid gives us a set of line in- 
tensities as a function of model parameters p = 
(Ni2 CO /dv,n(H 2 ),T k in,$A,X mo i/Xi2 CO ), which we 
then compare to our measured intensities x. To com- 
pare to column density per unit linewidth, we divide the 
measured intensities by 180 kms -1 , so that they are also 
per unit linewidth. The optical depth, and in turn the 
RADEX results, de pend only on the ra tio of column den- 
sity to line width. Wa rd et all (|2003l ) found linewidths 
of 180 kms -1 for the NE component and 160 kms -1 for 
the SW by resolving the structure in position- velocity di- 
agrams for their study from CO J = 1 ^ to J = 7 — s- 6, 
but we do not resolve the difference between the two and 
so we use the larger value. The Bayesian likelihood of 
the model parameters given the measurements is 



P(p\x) = 



P(p)P(x\p) 
P(x) 



(1) 



where P(p) is the prior probability of the model param- 
eters (see Section \3.2\i , P(x) is for normalization, and 
P(x\p) is the probability of obtaining the observed data 
set given that the source follows the model described by 
p. P(x\p) is the product of Gaussian distributions in 
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each observation, 



P(x\p) = TT y= — g exp 



2a 



(2) 



where oi is the standard deviation of the observational 
measurement for transition i and h(p) is the RADEX- 
predicted line intensity for that transition and model. 
For the total uncertainty, we take the statistical uncer- 
tainty in the total integrated intensity from the line fit- 
ting procedure and add 20%/10% calibration error for 
SSW/SLW in quadrature. To hnd the likelihood distri- 
bution of one parameter out of all of p, we integrate over 
all other parameters to find, for example, P(T&j„). 

3.1. Separate Components 

We divide the lines flu xes into two compone nts, one 
warmer and one cooler. iPanuzzo et al.l (|2010f) already 
showed that the high- J lines of 12 CO trace a warmer 
component than those transitions available from the 
ground. However, some of the mid- J lines (especially CO 
J = 4 — > 3) may have significant contributions from both 
components. To separate the 12 CO fluxes from each com- 
ponent, we follow an iterative pr ocedure. We first m odel 
the lowest three 12 CO lines from lWard et al.l (|2003D ; we 
take the sum of their measurements for the two observed 
lobes and scale the result by the ratio of their beam area 
to our 43" beam. The best-fit SLED is then subtracted 
from our SPIRE measurements, producing the black tri- 
angles in Figures [5] and [5] These triangles comprise the 
"warm component." The best fit warm SLED is then 
subtracted from the low-J lines, producing the asterisks 
in the aforementioned figure. These fluxes are refit to 
produce our results for the "cool component." 

We present a two-component model using just 12 CO, 
as well as one including our high- J lines of 13 CO along 
with the warm compone nt. We did not mo del the low- 
J 13 CO lines reported in iWard et al.l ([20031 ) due to the 
uncertainties presented in their Table 2 footnotes. We 
instead predicted the 13 CO spectrum of the cool compo- 
nent, given its best fit results and a 12 CO/ 1 3 CO ratio of 
of 35 (inbetween previously found 30 to 40 Wa rd et al.l 
2003), and found a very small contribution to the higher- 
J lines we are modeling. These small contributions are 
subtracted from the warm component, as with 12 CO in 
the previous paragraph. 

We also sought to include [C I], which was assumed 
to be associated with the cool component, due to 
the low excitation temperature (~ 30 K) derived from 
the line ratios (n u /ni — g u /giexp(—AT/T ex ), ni cx 
h [W / cm 2 } /(A^i)). Here, we use AT = 38.84, g u = 5, 
gi = 3, A u = 2.65 x 10 6 s -1 , A l = 7.88 x 10 8 s _1 . See Ta- 
ble [T] for the frequencies and unit conversion. However, 
this model produced some unphysical situations. We dis- 
cuss our findings and implications of them in Section [4.1l 

This analysis necessarily assumes that all of the line 
emission for a given component is coming from one por- 
tion of gas described in bulk by the model parameters. In 
reality, there is likely a variety of physical conditions, ex- 
isting in a continuum of parameters. However, the high- 
J SPIRE data does not provide justification for model- 
ing more than one warm component because the SLEDs 
are well-described by one component. We did attempt 



TABLE 3 
Likelihood Parameters Used 



Parameter 


Value a 


Units 


Line width b 


180 


[km s" 1 ] 


Abundance (Xi2qq / Xh 2 ) 


3.0 xlO" 4 


Angular Size Scale c 


17 


[pc/"] 


Emission Size 


43.0 


["] 


Length Limit 


900 


[PC] 


Dynamical Mass Limit 


2 xlO 9 


[M ] 



a Citations for parameters are in Section [32] 
b Used for scaling the line intensities. All other pa- 
rameters used for prior probabilities. 
c ^-region = T (Angular Size Scale [pc/"] X Source 
Size ["] / 2) 2 

a procedure to model multiple warm components of CO 
by first fitting the highest-J lines and subtracting the 
predicted line fluxes f or the mid-J lin e s. Su ch a proce- 
dure has been used in [Rangwal a et al.l ()2011[ ). for exam- 
ple. However, the predictions for the mid-J lines either 
matched or were an overestimate of the observed fluxes, 
leaving no second component to be modeled. A range of 
conditions is definitely present in the molecular gas, yet 
these two (warm and cool) components are dominating 
the emission, within the observational and modeling un- 
certainties. We note that, with regards to the different 
molecules/atoms being modeled here, all three species 
have similar profiles, as sho wn from the HIFI (h igher 
spectral resolution) spectra in lLoenen et al.l (|20100 . 

3.2. Prior Probabilities 

We use a binary prior probability, P(p), to indicate 
either a physically allowed scenario (P(p)=l) or an un- 
physical and thus not allowed scenario (P(p)=0). In 
other words, all combinations of parameters that are 
deemed physical based on the following three conditions 
were given equal prior probability, and all others are 
given zero prior probability. The conditions are as fol- 
lows: 

1. The total length of the column (L co i ) cannot exceed 
the length of the entire region. This assumes the 
length in the plane of the sky is the same as that 
orthogonal to the plane of the sky; we chose an 
upper-limit to the length of 900 pc because of the 
observed size (|Ward et al.1 120031. This is the most 
significant of all the priors, placing an upper limit 
on the column density and a lower limit on the 
density. This prior can be stated as 



Ai 



CO 



n(H 2 )V®AXi2 CO 
For the relative abundance X 



< 900 pc. 



(3) 



ror tne relative auunuaiice Ai2 p^ to molecular h y- 
drogen, we assume 3.0 xlO" 4 jWard et al.ll2003t ). 

2. The total mass in the emission region (M reg i on ) 
cannot exceed the dynamical mass of the galaxy. 
We use the expression 



M, 



region 



A r egionNi2 CO $ A l.5m H , 



Xl2 CO 



(4) 
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to calculate the mass in the emission region, where 
the 1.5 accounts for helium and other heavy ele- 
ments, and $a is the filling factor . We estimate the 
dynamical mass to be 2x 10 9 M (jWard et al.ll2003t 
iNaylor et al.1 12010). calculated using rotational ve- 
locity and radius. The other assumed parameters 
in the above expression are listed in Table [3] 

3. The optical depth of a line must be less than 100, as 
recommended by the RADEX documentation. The 
cloud excitation temperature can become too de- 
pendent on optical depth at high column densities, 
and so very high optical depths can lead to unreli- 
able temperatures. We found that in the presence 
of the other priors, this limit does not affect the 
likelihood results. 

3.3. Likelihood Analysis of the Map 

We run each map pixel through the aforementioned 
likelihood analysis independently. (Note that due to the 
beam size being larger than the pixels themselves, each 
pixel's data are not independent of its neighbors). We 
only model 12 CO in the spectral maps because they are 
of lower integration time and 13 CO cannot be reliably 
measured. We also cannot account for cool emission at 
different locations in our map. 

To be run through the likelihood analysis, a pixel was 
required to have both an SLW and SSW spectrum and at 
least 5 12 CO lines with S/N > 10 (the convolution tends 
to increase the signal/noise of each pixel). 90 pixels met 
this requirement (98 pixels would meet the requirement 
of 5 12 CO lines with S/N > 3, so little would be gained by 
going to lower S/N). Additionally, we do not find statisti- 
cally significantly different results requiring only 3 lines, 
the minimum with which we could reasonably model the 
emission; to some extent, this requires at minimum the 
J = 6 — > 5 transition, which means we will be tracing 
higher temperatures. 

4. MODELING RESULTS AND DISCUSSION 
4.1. Physical Conditions: Deep Spectrum 

We present two different versions of our likelihood anal- 
ysis for the deep spectrum: one using only 12 CO, and one 
using 12 CO and 13 CO ("multiple molecule"). The moti- 
vation behind this is to investigate two questions: does 
the addition of different species change the modeled pa- 
rameters of the gas, and/or does it better constrain the 
parameters? Both versions contain a warm and cool com- 
ponent. The modeling assumes all of the emission in a 
given component is coming from the same homogeneous 
gas, and by comparing these models we will investigate 
the validity of this assumption in this subsection. 

The results for each of these versions are presented in 
Tables 0] and O respectively. Figures [5] and [5] show the 
input SLED as well as the best-fit model results. The 
primary results (temperature, density, column density, 
and filling factor) are displayed graphically in Figures [7] 
and [8] Secondary parameters, which are calculated from 
the aforementioned primary results, are displayed in Fig- 
ures [9] and [TU] These include the pressure (the product of 
temperature and density) and the beam-averaged column 
density ((-/Vi2 CO )j the product of column density and fill- 
ing factor). We note that in the results, the parameters 
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Fig. 5. — Bayesian Likelihood Analysis, Spectral Line Energy 
Distributions, 1 CO Only. Asterisks represent the cool component, 
with its best fit SLED ("4D Max" column in Table gjl shown by 
a dotted line. Triangles represent the warm component, with its 
best fit SLED shown by a dashed line. The total measurements are 
shown with diamonds with their associated error bars. The total 
of both components is the solid line. 
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Fig. 6. — Bayesian Likelihood Analysis, Spectral Line Energy 
Distributions, including 13 CO. Each color is a separate species: 
black for 12 CO, red for 13 CO (only warm component). The total 
fluxes are shown by diamonds, but their separate components are in 
asterisks/triangles for the cool/warm components. Best fit SLEDs 
("4D Max" column in Tabled are shown with dotted/dashed lines 
for cool/warm components. The total SLED, shown with a solid 
line, is the sum of the two components. 
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TABLE 4 
Likelihood Results: 12 CO Only 





Integrated Likelihood: Cool 


Cool 


integrated Likelihood: Warm 


Warm 






Median 


1 Sigma Range 


ID Max 


4D Max 


Median 


1 Sigma Range 


ID Max 


4D Max 




Tkin 


40 


12 - 472 


13 


63 


414 


335 - 518 


447 


447 


[Kj 


n(Ha) 


10 3.23 


1Q2.36 _ 1Q 4.81 


10 2.20 


10 3.40 


10 3.98 


10 3.53 _ 10 4.21 


10 410 


1Q 4.10 


[cm 3 ] 


N co 


1Q19.03 


IO 18 - 38 _ io 19 - 56 




■LQ18.56 


10 18.12 


10 1706 — io 1909 


10 18 -56 


10 17.96 


[cm -2 ] 


*A 


10 -0.92 


10 -1.22 _ 10 -0.57 


10 -0.97 


10-0-9° 


1Q -1.50 


10-i- 9 6 - IO" 52 


10 -1.88 


1Q -1.28 




P 


10 5.11 


10 4.55 _ 10 S.76 


10 5.24 


10 5.24 


10 6.61 


1 8-18 _ 1 6.80 


10 6.75 


10 6.75 


[K cm" 2 ] 


< Nco > 


10 18.15 


10 17 - 63 — io 18 - 84 


10 17.78 


10 17.78 


10 16.67 


10 16 44 - io 17 - 20 


-LQ16.58 


10 16 - 58 


[cm -2 ] 


Mass 


10 7.67 


107-16 _ 10 8.36 


10 7.31 


10 7.31 


10 6.20 


1 5- 97 _ 1 6 72 


10 6.11 


10 6.11 


[M ] 



Note. — ID Max refers to the maximum likelihood of that parameter after integrated over all the other parameters (the 
mode of the likelihood distributions). 4D Max refers to the single most probable grid point (mode of the entire multi-dimensional 
distribution). Median, 1 sigma lower and upper values refer to the integrated distribution, when the cumulative distribution 
function equals 0.5, ~ 0.16 and ~ 0.84, respectively. Note that the lD Max may be outside the 1 sigma range because of 
asymmetry in the integrated likelihood. P and (Ni2qq) are calculated from the 2D distribution of and n(H2) (or <I>a and 
Ni2qq) , which is why we give the same value in both the ID Max and 4D Max columns. 



TABLE 5 

Likelihood Results: 12 CO and 13 CO 





Integrated Likelihood: Cool 


Cool 


Integrated Likelihood: Warm 


Warm 






Median 


1 Sigma Range 


ID Max 


4D Max 


Median 


1 Sigma Range 


ID Max 


4D Max 






35 


12 - 385 


14 


158 


436 


344 - 548 


447 


501 


[K] 


n(H 2 ) 


10 3.44 


10 2.48 _ 10 5.15 


10 3.oo 


10 3.10 


10 3.58 


IO 3 " - io 3 - 96 


10 380 


10 3.40 


[cm" 3 ] 


N co 


10 19.25 


10 18 -69 _ jo 19 - 70 


1 0l 93 6 


-Lq18.96 


1Q 19.02 


10 1819 — io 19 - 51 


10 19.16 


10 19.26 


[cm -2 ] 


$ A 


10-1-19 


10-1.82 _ 10 -0.83 


1Q-1.28 


10 -1.3S 


lQ-1.88 


10-2 06 _ 10-1-50 


IO" 1 ' 88 


10 -1.88 




P 


10 5.24 


10 4 73 - io 6 - 18 


10 5.24 


10 5.24 


10 6.23 


1Q5.S0 _ io 6 - 60 


10 6.61 


10 6.61 


[K cm- 2 ] 


< N co > 


1Q,18.09 


10 17.59 _ 1Q 18.71 


10 17.78 


1Q 17.78 


10 17.12 


1Q16.69 _ io 17 - 59 


10 16.76 


10 16.76 


[cm -2 ] 


Mass 


1Q 7.62 


10 7 11 - io 8 - 23 


1Q 7.31 


10 7.31 


1Q 6.65 


1Q 6.22 _ 1Q 7.12 


10 6.28 


1Q 6.28 


[M ] 


-^"l3co/ N C o 










10-1.58 


10-i-TU _ i -i- 4 e 


10 -1.SO 


10 -i.so 




N 13co 










10 17.44 


IO 16 62 _ io 17 - 92 


10 17.62 


10 17.62 
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Fig. 7. — Bayesian Likelihood Analysis, Primary Parameter Re- 
sults, 12 CO only. Each color represents a separate component; blue 
for cool, red for warm (see Section[3jl. Dashed/dotted vertical lines 
indicate the median/4D maximum of the distribution (see TableQ. 
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Fig. 8. — Bayesian Likelihood Analysis, Primary Parameter Re- 
sults, including 13 CO. Each color represents a separate component; 
blue for cool, red for warm (see Section [2J- Dashed/dotted verti- 
cal lines indicate the median/4D maximum of the distribution (see 
Table Q. 



of the most likely grid point ( "4D Max" ) is not necessar- 
ily the same as the median or the mode ("ID Max") of 
the integrated likelihood distributions. The "4D Max" is 
describing one specific point, but the median, ID Max, 



and associated error range are representative of the larger 
likelihood across all other parameters in the grid. 

We first compare our two models, 12 CO only vs. mul- 
tiple molecule. These are Tables|4]vs. [U Figures[7]vs. [51 
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Fig. 9. — Bayesian Likelihood Analysis, Secondary Parameter 
Results, 12 CO only. Top: Each color represents a separate com- 
ponent; blue for cool, red for warm (see Section [3j Dashed/dotted 
vertical lines indicate the median / 4D maximum of the distribution 
(see Table |4j. Bottom: 2D distributions for the pairs of primary 
parameters from which the secondary parameters (top) were de- 
rived; colors correspond to component. Diagonal lines indicate 
values of the top parameters (pressure on left and beam-averaged 
column density on right). Contour levels are 20, 40, 60, and 80% 
of the maximum likelihood. 

and Figures [9] vs. QUI 

The addition of 13 CO to the warm component does not 
significantly change the temperature, but it does increase 
the likelihood of lower densities. It also decreases the 
likelihood of lower column densities. The consequences 
of these changes can be seen in the pressure and mass dis- 
tributions in Figures l9l and ITUl Adding 13 CO increased 
the likelihoods of the "shoulders" of these distributions; 
the lower half of pressure, and the upper half of mass. 
An examination of the contour plots in the bottom half 
of these figures illustrates why, statistically. In the 12 CO 
only model, though column density and filling factor are 
not well constrained independently, they are highly cor- 
related; their contours run along an almost constant line 
of beam-averaged column density. Adding 13 CO intro- 
duced the Xi3co/Xi2co parameter, which also impacts 
the absolute flux level of the models, like column den- 
sity and filling factor. The result are likelihoods that are 
more constrained but not as highly correlated with one 
another. Therefore the mass distribution is wider. In 
the 12 CO-only model, the mass of the warm component 
is about 3.4 /6.3% (median/4D Max) the mass of the cool 
component. iRigopoulou et al.l (|2002f) noted that "warm" 
gas is generally around 1 to 10% of total gas mass for 
starburst galaxies, so this is about as expected. The ad- 
dition of 13 CO creates somewhat overlapping likelihood 
distributions for mass (the 1-sigma ranges are just touch- 
ing), but the median and 4D Max now warm/cool ratios 
of 11% and 9%, respectively. One factor that may con- 
tribute to wider distributions is overestimated error; the 
error bars are dominated by our 20% calibration error, 
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Fig. 10. — Bayesian Likelihood Analysis, Secondary Parame- 
ter Results, including 13 CO. Top: Each color represents a sep- 
arate component; blue for cool, red for warm (see Section [3} 
Dashed/dotted vertical lines indicate the median/4D maximum of 
the distribution (see Table [4}. Bottom: 2D distributions for the 
pairs of primary parameters from which the secondary parameters 
(top) were derived; colors correspond to component. Diagonal lines 
indicate values of the top parameters (pressure on left and beam- 
averaged column density on right). Contour levels are 20, 40, 60, 
and 80% of the maximum likelihood. 

not statistical error. 

After this point , we compare frequently to 
IPanuzzo et al.l (|2010f) . The major factor responsi- 
ble for the differences, just modeling 12 CO alone, is 
the shape of the CO SLED at those lines with upper 
rotational number greater than J=8. We also explicitely 
subtracted the cool co mponent's cont r ibutio n from 
12 CO J = 4 -)• 3, whereas IPanuzzo et~aT1 (|2010D simply 
unde rpredicted the total flux. We will also compare 
with ILoenen et al.l (|2010D , another high- J CO study of 
M82, in Section E31 

Our results are similar to IPanuzzo et al.1 ((2010) , who 
found that these high-J CO lines trace a very warm gas 
component that is separate from the cold molecular gas 
traced by those lines below J — 4 — > 3. Our best-fit 
temperature of the 12 CO only model (at 447 K, Ta- 
ble HJ is close to their value at 545 K, but the overall 
likelihood for temperature, integrated over all parame- 
ters, yields a slightly lower 414 K. Given the size of the 
uncertainty (335-518 K) in the parameter, the two dis- 
tributions are very similar, and therefore the result is 
not significantly different. Such warm gas has also been 
traced in the S(l) and S(2) transition s of H 2 , at 450 K 
with the Infrared Space Observatory (jRigopoulou et al.l 
|2002|) and 536 K w ith the Spitzer Infrared Spectrograph 
(|Beirao et aHl2008ft . 

We do find a slightly higher density than lPanuzzo et all 
(2010), with our best-fit value of 10 41 compared to 
10 37 cm -3 , though the integrated likelihood distribu- 
tions do overlap (see Figure [7]). However, the temper- 
ature and density are degenerate; higher temperatures 
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and lower densities may produce the same fluxes as lower 
temperatures and higher densities. Their product, the 
pressure, is better constrained. We seem to have col- 
lapsed/constrained the p ressure distribu t ion to the upper 
half of that presented in iPanuzzo et al.1 (|2010( ) . 

The co lumn density is not a s well constrained as pre- 
sented in IPanuzzo et al.l (|2010() ; we found that they had 
an error in calculating the expected fluxes of the higher- J 
lines for lower column density values. We have recalcu- 
lated the fluxes for those column densities, and we find 
that in fact when properly calculated these column densi- 
ties have a non-zero likelihood. In the 12 CO only model, 
the column density itself is not constrained. However, 
the column density and filling factor are degenerate, so it 
is their product (beam- averaged column density, (N co )) 
that is better constrained. Our best-fit value is 10 16 6 
cm~ 2 . 

The total mass in the beam can be calculated using 
Equation 2] (and is presented as the top y-axis in Fig- 
ures M and QUI upper right). As previously discussed, the 
12 CO only model produces the expected result of less gas 
mass in the warm component. In the cool component we 
find a best-fit mass of 2.0 xlO 7 M© (median 4.7 xlO 7 
M Q ). This is small er than the 2.0 xl O 8 M Q traced by 
the LVG analysis of IWard et all ([200l with lower-J CO 
lines. This difference is due to the fact that we subtract 
the contribution to the low-J flux from the warm com- 
ponent; in our initial modeling of the cold component, 
before this subtraction, our best-fit mass is 9.8 x 10 7 M© 
(with a range from 0.2 to 2.2 x 10 8 M ). The warm com- 
ponent is a smaller fraction of the gas, with a best-fit of 
1.3 xlO 6 Mq, about 6.3% the mass of the cool compo- 
nent. 

The 12 CO/ 13 CO relative abundance is also a free pa- 
rameter in our multi-species model; we find a best-fit 
value of about 32, similar to the 40 and 30 found previ- 
ously for the NE and SW lobes, respectively ()Ward et al.l 
l2003h . 

As mentioned in Section [3. 11 we also attempted to in- 
clude our two [C I] lines with the cool component. We 
do not present the tables for this model because the 
mass distributions of the warm and cool components be- 
came overlapping, indicating the same amount of mass 
in both components, an unphysical situation. Addition- 
ally, the derived relative abundance of [C I] to 12 CO 
was unusually high. We found a ratio of 0.48 to 3.3, 
which i s higher than iWhite et al.l (|1994 average value 
~ 0.5 ), ISchilke et all ((199J, 0.1-0.3), and IStutzki et al] 
((19971 0.5) using other methods. Before subtracting the 
warm component's contribution to the 12 CO flux (when 
we just fit all of the low-J 12 CO flux and [C I]), we find ra- 
tios more consistent with these values (best-fit 0.4, range 
0.09 to 1.23). These two problems could be indicating 
that the assumption of CO and [C I] coming from the 
same component is flawed. The column density, temper- 
ature, and mass developed somewhat of a double-peaked 
structure; specifically, the addition of [C I] increases the 
likelihood of lower column densities and masses, but does 
not eliminiate the previous likelihood peak. 

It is unclear how much of the molecular CO and atomic 
C are truly cospatial and therefore how physical our re- 
sults for modeling them all together as o ne bu lk gas com- 
ponent may be. IPapadopoulos et all (|2004D presented 



results which argue that [C I] and CO are cospatial and 
trace the same hydrogen gas mass ([C I] doing so better 
than CO in many conditions). This conflicts with the 
theoretical picture of gas clouds (especially in PDRs) as 
a structured transition between molecular, atomic, and 
ionized gas, but new observational and theoretical evi- 
dence indicates t he types of gas are not so distinct (see 
references within IPapadopoulos et all 2004 ). For exam- 
ple, iHowe et al.l ([2000D and iLi et al.l (|2004l) have found 
[C I] to trace 13 CO well. If the ISM is clumpy, well 
mixed, and dynamic, the [C I] and CO may be cospa- 
tial averaged over large scales. Strong stellar winds (and 
possibly the interaction with M81) could be contributing 
to the dynamic nature of the gas, so it is not unreason- 
able to believe that the gas has not achieved the simple 
layered pattern. Though the SPIRE FTS cannot resolve 
the two separate velocity components of M82, HIFI can, 
and observations of these two [C I] lines indicate a gen- 
erally similar shape to the CO lines, namely two Gaus- 
sians with the SW c omponent demonstrating higher flux 
(ILoenen et alj|2010l). 

iWolfire et all (|2010f l presented a model PDR which 
shows a cloud layer traced by atomic (not-yet-ionized) 
[C I] where the hydrogen is still molecular due to self- 
shielding effects. This "dark molecular gas" (called so 
because it is not traced by CO) would be less-shielded 
and at a warmer temperature than the inner-most cloud 
layer of CO. It is possible that [C I] and CO are some- 
what cospatial yet somewhat segregated as in the "dark 
molecular gas model." Our analysis is consistent with a 
picture in which the [C I] and 12 CO do not completely 
overlap spatially. 

The [C I] J = 1 — > emission can also be used 
to est imate the total hyd r ogen mass using Equation 
12 of Papadopoul os et all (|2004l ). Using the median 
Xf CI] /X H2 of 1.5 xl0~ 4 from Xjc i]/Xi 2co = 0.5 
([White et alj|1994t IStutzki et alj|1997l ) and our assumed 

Xi2 CO /Xh 2 , we find a total gas mass of 4.4 x lO 7 ^^ 1 
M . Qio is the ratio of the column of the J = 
1 — s- emis sion to the total [C I] c olumn (see Ap- 
pendix A of IPapadopoulos et all 12004) . which depends 
on the excitation conditions of the gas; for Qio ~ 0.5 
(|Papadopoulos k. Grevell2004f ). the gas mass is 8.8 x 10 7 
Mq but is uncertain by modeled uncertainty in X\ct\ 
alone. This method of estimating the mass is higher than 
total mass estimate of the cool component described ear- 
lier in this section. [C I] may be coming from a range 
of temperatures, but with only two lines, we cannot sort 
that out. 

4.2. Physical Conditions: Map 

Results for all of the same parameters for the deep 
spectrum were produced for each of the pixels in the 
map (that met the criteria in Section [373]). We find that 
the beam size of SPIRE cannot resolve the structure in 
M82 as has been done with interferometric maps or high 
spectral resolution observations (which can resolve the 
velocity components of the NE and SW lobes). Because 
of the degeneracy between temperature/density and col- 
umn density/filling factor, we present results of their 
products, pressure and beam-averaged column density 
in Figure [TT1 (Ni2 CO ) shows a radially decreasing trend, 
roughly corresponding to the decrease in the beam pro- 
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file (plotted with a solid black line), implying an obser- 
vational effect due to source-beam coupling. We would 
expect an off-centered beam to be able to probe the pres- 
sure of the central region (because the relative ratios of 
the SLED would be preserved), so the lack of a radial 
trend also implies that we are not measuring different 
areas of M82 in each map pixel. This indicates that 
the SPIRE FTS map cannot resolve M82's structure, 
and therefore the single "deep" spectrum is an adequate 
representation of the galaxy as a whole. Our results in 
Section FOl are descriptive of the bulk properties of the 
galaxy and we do not see trends on the scale of our map. 
This analysis is separate from the dust temperature gra- 
dients (which we also see in the co ntinuum gradients o f 
our map spectra) found in M82 by iRoussel et al.l (|2010f ) 
and indicates that the dust and molecular gas are not 
coupled. 

One difficulty with the map is it has lower signal/noise; 
it must also be convolved up the largest beam size (43", 
like the deep spectrum). However, as mentioned in Sec- 
tion 14.11 it is the highest- J fluxes in the SSW detector 
that largely constrain the results of the likelihood for the 
deep spectrum. Therefore, we also attempted to model 
the map with just those lines with upper-J level of 9 or 



higher, without convolving. These lines all were mea- 
sured with a beam FWHM of ~ 19", offering higher res- 
olution. However, these lines are also weaker, and with 
fewer, lower signal/noise lines this method does not con- 
strain any parameters as well. 

Though the off-axis pixels may not provide new infor- 
mation about the physical trends of the galaxy, we can 
compare the deep spectrum to the center pixel of our 
map as a test of the source-beam coupling corrections 
(see Section 12. 2|) . The results of the cent ral p ixel are 
very similar to those presented in Section 14.11 though 
the integrated likelihoods are generally wider (the pa- 
rameters are less constrained). This is partially, but not 
entirely, due to larger error bars on the SSW lines. 

4.3. Molecular Line Absorption 

In addition to the Bayesian likelihood modeling, we 
can also briefly discuss the absorption lines presented in 
Table ffl 

4.3.1. Hydrogen Fluoride (HF) 

Hydrogen fluoride (HF) is potentially a sensitive probe 
of total molecular gas column, because the HF/H2 ra- 
tio is more reliably constant than 12 CO/H2 because the 
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for mation of HF is do minated by a reaction of F with 
H 2 (jMonie et al.l 120111 ). Furthermore, HF J = 1 ->• 
is generally seen in absorp tion because of its high A- 
coefhcient, 2.42 xlCT 2 s _1 (|Monje et al.ll2011[ ). Assum- 
ing all HF molecules are in the gro und state (generall y 
true in the diffuse and dense ISM, iMonje et al.l l20lTI ). 
the HF J = 1 — > line yields the optical depth simply as 
r = —ln(Fi/F c ), where Fi/F c is the line-to-continuum 
ratio. In the case of HF, we mask out a nearby water 
emission line (1226 to 1229 GHz), though because the 
signal in each wavenumber bin is not independent due to 
the ringing, this can introduce added uncertainty. There- 
fore the following discussion is meant to be approximate. 
We estimate the HF column density using 



8ngi 



(5) 



where g u = 3 and gi = 1 . This implies J rdv = 
4.16 x W- 13 N(HF)cm 2 km s" 1 . The HF line occurs 
in a part of our spectrum with some noticeable struc- 
ture in the continuum (see Figure |3l third panel, around 
1250 GHz). If we only integrate the 6 GHz surrounding 
the line, we find a column density of HF of 6.61 xlO 13 
cm" 2 . Expanding the range over which we integrate in- 
creases the derived column density, but this may be due 
to other features in the spectrum, and so we consider 
our derived value a lower limit. Assuming a pred icted 
abundance of HF of 3.6 x 10" 8 (jMonie et al.ll201lD . this 
corresponds to a molecular hydrogen column density of 
1.84 x 10 21 cm -2 . The column density derived from this 
line is similar to that of (N) of the cool component of 
12 CO. However, there are still some uncertainties to this 
calculation. We are only able to see the HF in front of 
the continuum emission, and therefore we are not prob- 
ing the total column density. Higher spectral resolution 
could reveal the extent of spatial colocation of HF with 
CO. There are also uncertainties associated with either 
molecular abundance assumed and whether or not all HF 
molecules are truly in the ground state. 

4.3.2. Water and Water Ion (H 2 {+) ) 

Water is fundamental to the energy balance of col- 
lapsing clouds and the subsequent formation of stars, 
planets, and life. Many Herschel key programs are cur- 
rently studying water and chemically related molecular 
species in a variety of conditions. An excellent summary 
of water chemistry i n star forming regions is available in 
Ivan Dishoeck et all (|2011l ). 

IWeifi et all (|2010l ) studied the low-level water transi- 
tions in M82, and they detect the ground-state 0-H2O 
emission in two clearly resolved components, which we 
do not. With the 41" beam, the two components add 
to 370 ± 44 Jy km/s beam" 1 , well below our thresh- 
old of detection, as can be seen by examining Table [TJ 
Though we do not detect that line, we do det ect four new 
water lines in addition to those presented in IWeifi et"aTI 
(2010): two p-H 2 (752 and 1229 GHz) and two o-H 2 
(1097 and 1153 GHz). Combined with their ground-state 
transition of 0-H2O, we add to the picture of the water 
excitation in M82. 

Our g round-state l ines i ndicate similar column den- 
sities as IWeifi et all ((2010), within a factor of 2. Us- 



ing Equation [5] (low-excitation approximation), we find 
column densiti es of p-H 2 Q and o-H20 + of ~ 4 x 10 13 
cm" 2 , whereas IWeifi et all (|2010f ) finds 9.0 and 2.2 xlO 13 
cm" 2 , respectively. They found that the water absorp- 
tion comes from a region northeast of the central CO 
peak; shocks related to the bar structure of M82 could 
be releasing water into the gas phase at such a location. 
The fact that the water comes from a lower column den- 
sity region seems to contradict the existence of a PDR, 
which would require high column densities to shield water 
from UV dissociation; however, the relative strength and 
similarity of absorption profiles of o-H20 + compared to 
P-H2O indicates some ionizing photons (see their work 
for complete interpretation). Though these transitions 
are tracing a different region than CO, they add to the 
picture of a complicated mix of energy sources present 
in the g addressed in Section 14.41 Models of wa- 

ter emis sion from shocks have been invest igated by oth- 
ers (i.e. iFlower fc Pineau des Foretsll2010l ). but detailed 
modeling of the water spectrum of M82 is outside the 
scope of this work. 

4.4. Gas Excitation 

At the high temperature of the war m component, the 
coolin g will be dominated by hydrogen. iLe Bourlot et al.l 
(1999) modeled the cooling rates for H2, and made their 
tabulated rates available with an interpolation routine 
for desired values of density, temperature, ortho- to para- 
H2 ratio, and H to H2 density raticFl. For our best-fit 
temperature a nd density, assuming n( H) /n(H2 ) = 1 (rec- 
ommended by ILe Bourlot et al.l 119991 for PDRs), this 
corresponds to a cooling rate of 10 -22 - 54 erg s" 1 per 
molecule, or 3 L Q /M Q (using o/p H 2 = 1, though the 
number is only ~ 3% lower for o/p H2 = 3). Given 
the warm mass of 1.3 x 10 6 M© (using the 12 CO model 
for the rest of this section) , that implies a hydrogen lu- 
minosity of 3.9 xlO 6 Lq. The total observed hydrogen 
luminosity thus far has been higher t han this number; 
addin g the luminosities presented in iRigopoulou et al.l 
(2002J) and correcting for extinction as in IDraind (11989). 
we find a total of 1.2 xlO 7 L Q in the (0-0)S(0)-S(3), S(5), 
S(7), and (1-0)Q(3) lines. However, some of these hy- 
drogen lines are tracing lower or higher temperature gas. 
We note that the mass range within one standard devi- 
ation of our likelihood results for the warm component 
is 0.93-5.2 xlO 6 M Q , which corresponds to a predicted 
luminosity of 0.28-1.6 xlO 7 L©, encompassing the mea- 
sured hydrogen luminosity. 

There are a few possibilities for the source of the exci- 
tation of the gas: X-ray photons, cosmic rays, UV exci- 
tation of PDRs and shocks/collisional excitation. Hard 
X-rays from an AGN have already been ruled out by 
others in the literature due to the lack of evidence for 
a strong AGN and low X-ra y luminosity (1.1 xlO 6 L©, 
IStrickland fc Heckmanll2007l) . 

The CO emission from M8 2 has previously b een in- 
terpreted using PDR models. iBeirao et al.l (|2008f ) noted 
with the Spitzer Infrared Spectrograph that the H2 emis- 
sion is correlated with PAH emission, indic ating that it is 
mainly excited by UV radiation in PDRs. lLoenen et al.l 
(2010) combined HIFI data with ground based detec- 

12 http:/ /ccp7. dur.ac.uk/cooling_by_h2/ 
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tions in order to model 12 CO J = 1 to J = 13 12 
and 13 CO J = 1^0toJ = 8^7. They reproduced the 
measured SLEDs with one low-density (log(n(-ff2))=3) 
and two high-density (log(n(-ff2))=5,6) components with 
relative proportions of 70%, 29%, and 1%, respectively. 
The low-density component is largely responsible for the 
low-J emission, while the highest-density component is 
responsible for the highest-J emission. These high densi- 
ties are not consistent with the results of our likelihood 
analysis detailed in Section 14. 1( the likelihood of solu- 
tions for the warm compo nent at log(n(ff2))>5 is essen- 
tially zero. We note that lLoenen et al] (|2010l )'s Figure 
3 sho ws the consiste n cy be tween HIFI and SPIRE fluxes 
from iPanuzzo et al.l (|2010f) . In other words, the differ- 
ence is not due to discrepant line fluxes, but different 
models (PDR vs. CO likelihood analysis). 

There are two major differen ces in the approach of this 
work and lLoenen et all ()2010h . First, the order in which 
we approach the problem is different. We first analyze 
the CO excitation using likelihood analysis to determine 
the physical conditions. Once we have these conditions, 
we then look to the possible energy sources based on the 
conditions wc have already derived, instead of first see- 
ing under which conditions a certain energy source fits. 
Second, we look beyond the one best fit solution: in ad- 
dition to presenting the best-fit solution, our likelihoods 
analyze the relative probabilities for a larger parameter 
space. 

We also attempted to reproduce our deep SLED with 
various PDR models. iMeijerink et al.l (j2006[ ) have added 
to their PDR and XDR models to include enhanced cos- 
mic rays (at a rate of 5 xlO" 15 s" 1 ), near our assumed 
rate discussed later in this section. Such models are cur- 
rently available for incident flux of log(Go) of 2-4 (Go = 
1.6 xlCT 3 erg cm" 2 s" 1 ) for log(rc(if 2 ))=3 and log(G ) 
of 3-5 for log(n(i?2)) of 4 and 5. By examining the ra- 
tios of 12 CO J = 9 — » 8 with all higher-J lines (those 
largely driving the likelihood results, and also measured 
from similar beam sizes), none of the available 9 PDR 
models are an ideal match, but the PDR scenario for 
\og(n(H 2 ))=3 7 log(Go)=3 is the best match. Figure 1121 
compares the predicted and observed ratios. The ra- 
tios used are without source-beam coupling correction 
because the SSW already has similar beam sizes, but 
the ratios with beam correction are within 4-8% of those 
presented in Figure [12] 

Additionally, we used a higher-resolution (in den- 
sity and incident flux) grid of 12 CO PDR models 
(jWolfire et al.l l201?fl . The same line ratios previously 
discussed (those shown in Figure [12]) can only be repro- 
duced by higher densities. For ( J = 9 — > 8)/ (J = 10 -> 9), 
the observed ratio is only found for log(n.(i?2)) > 4.5, 
log(Go) > 2, and by (J = 9 -> 8)/(J = 13 -> 12), only for 
log(n(JT 2 )) > 5, log(Go) > 2.5. 

To summarize, current PDR models can only explain 
the observed high- J 12 CO emission with densities higher 
than those indicated by the likelihood analysis (even 
when all priors are excluded). The cosmic-ray enhanced 
PDR models, though sparse, can come closer to repro- 
ducing the high-J line ratios at a lower density, though 
these are also below our likelihoods. There is also ev- 
idence that shocks enh ance high-J lines far more than 
PDRs (TPon et all [2012(1 . In their models, almost all of 
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the emission from the lowest- J lines came from unshocked 
gas, and most of the emission above J = 7 — > 6 was 
from shocked gas. The combination of shocks and PDRs 
could be responsible for the high-J CO lines observed, 
while PDRs alone are adequate to explain the lower-J 
CO lines and PAH emission. In summary, these high-J 
lines are not consistent with current PDR models, but 
improved modeling of shocks and the effects of cosmic 
rays on PDRs may help explain their emission. 

Cosmic rays (CRs) are another possibility for the ex- 
citation. The Very Energetic Radiation Imaging Tele- 
scope Array System (VERITAS) Collaboration recently 
reported that the cosmic-ray density in M82 is about 
500 times the average Galactic (M ilky Way) density 
(| VERITAS Collaboration et al.ll2009h . By using the cos- 
mic ray ionization rate of the Galaxy (2-7 xl0~ 17 s" 1 , 
iGoldsmith fc Langer|[l978t Ivan Dishoeck fc Black|[T986T) . 
multiplied by 500, and then multiplied by the average 
energy per ionization (20 eV), one finds an energy depo- 
sition rate of 2-7 xlO" 13 eV/s per H2 molecule in M82. 
This implies a heating rate of 0.03 to 0.12 L0/M0, less 
than 5% of the required molecular hydrogen cooling rate. 
Thus, cosmic rays alone cannot excite the molecular gas. 

Turbulent heating mechanisms may also play a role in 
M82. From lBradford et al.l (|2005[ ). the turbulent heating 
per unit mass can be expressed in L Q /M as 



1.10 



25 km s 



1 pc 

a7 
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where v rm s is the turbulent velocity and A^ is the typi- 
cal size scale of turbulent structures. We use the Jeans 
length for this size scale, calculated from the parame- 
ters of the likelihood results (Section 14.11) . which is 0.9 
pc. Given this size scale, the observed cooling rate could 
be replicated with a turbulent velocity of 33.7 kms" 1 . 
This would imply a velocity gradient of approximately 
37.5 kms" 1 pc" 1 . When we calculate the velocity gradi- 
ent using our model results (dv/dr = Av n{H2) / N {H2)) , 
we find a 68% confidence lower limit of 16 kms" 1 pc" 1 ] 
(4 kms~ 1 pc" 1 if we use the model with 13 CO) , but 
the upper limit is unphysically high. The velocity re- 
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quired for turbulent hea ting seems reasonable in context 
of our likelihood results. iPanuzzo et a l. ( 2010) used their 
calculated velocity gradient of 35 kms -1 pc _1 to deter- 
mine that they could match the heating required with 
a sizescale of 0.3 to 1.6 pc. These velocity gradients 
seem large compare d to Galactic star-forming sites (e.g. 
Ilmara fc Blitd [201 If ) . but M82 is known to have p ower- 
ful stellar winds. According to iBeirao et al.l (|2008l ). the 
starburst activity has decreased in the past few Myr, and 
this appears to be evidence of negative feedback (by stel- 
lar winds and supernovae), because M82 still has a large 
reservoir of gas available for star formation. Therefore, 
there is evidence for tu rbulent heating mechanism s being 
in place. Additionally, iDownes fc Solomonl (|1998| ) found 
high turbulent velocities (30-140 kms -1 ) in models of 
extreme star-forming galaxies. 

None of the possibilities described seem to provide 
enough heating by themselves, with the exception of 
turbulent heating, which is based on a few approxima- 
tions and assumptions. Likely, there is a combination 
of factors, namely PDRs and shocks/turbulent mecha- 
nisms. Such a situation has also been seen in other 
submillimeter-bright galaxies, discussed next. Interest- 
ingly, even a more quiescent galaxy like NGC 891 requires 
a combination of PDRs and sho cks to explain mid - J CO 
transitions (J = 6->5, J = 7->6. lNikola et al.ll201lD . 

4.5. Comparison to Other Starburst and Submillimeter 

Galaxies 

Because only the first few lines in the CO ladder are 
easily visible from the ground for nearby galaxies, the 
high- J lines detected by Herschel represent new territory. 
Therefore, while adequate diagnostics of high-J CO lines 
are still being developed, it is useful to compare to other 
submillimeter-bright galaxies. 

Mrk 231 contains a luminous (Seyfert 1) AGN. It also 
shows a strong high-J CO ladder, such that only the 
emission up to J — 8 — > 7 is explained by UV irradia- 
tion from star formation. Their high-J CO luminosity 
SLED however, is fiat (though ours for M82 are stronger 
than predicted, they are still decreasing with higher- J), 
van der Werf et all ([2010h can explain this trend with ei- 
ther an XDR or a dense PDR. An additional difference 
between M82 and Mrk 231 is that OH+ and H 2 0+ are 
both seen in strong emission in Mrk 231 (instead of ab- 
sorption), indicative of X-ray driven chemistry. Mrk 231 
is also more face-on than M82. 

The FTS spectrum of Arp 220 has many features not 
present in M82, such as strong HCN absorption, P-Cygni 
profile s of OH + , H2O+ and H2O, and evidence for an 
AGN (Rangwala e t al.l 120111 ). CO modeling, similar to 
the procedure done in this work, also indicates that the 
high-J lines trace a warmer component than low- J lines 
(~ 1350 K). Mechanical energy likely plays a large role in 
the heating of this merger galaxy as well. Though M82 
has an outflow, it is not detected in P-Cygni profiles of 
the aforementioned lines. 

The redshift of HLSW-01 (jConlev et alJl20Tl allows 
the CO J = 7— >6 to J = 10— >9 lines to be ob served from 
the g round, as has been done with Z-Spec (jScott et al.l 
l20Tll) . Unlike M82 (and others), the known CO SLED 
from J = 1 — > and up can be described by a single 
component at 227 K (1.2 x 10 3 cm~ 3 density). If the 



velocity gradient is not constrained to be greater than 
or equal to that corresponding to virialized motion, the 
best fit solution becomes 566 K (0.3 x 10 3 cm -3 density), 
closer to our temperature. We chose to exclude this prior 
due to uncertainties in the calculation of velocity gradi- 
ent related to M82's turbulent morphology. HLSW-01 
appears to be unique in that a cold gas component is 
not required to fit the lower- J lines of the SLED, though 
two-component models can find a best-fit solution with 
a cold component. 

In summary, M82 (like Arp 220 and HLSW-01) does 
not have the high CO excitation dominated by an AGN 
as seen in Mrk 231. Therefore, in addition to distinguish- 
ing between PDRs and shocks, high-J CO lines may also 
be used to indicate XDRs. 

5. CONCLUSIONS 

We have presented a multitude of molecular and atomic 
lines from M82 in the wavelength range (194-671 /zm) 
accessible by the Herschel-SPIRE FTS (Tabled)). After 
modeling 12 CO, 13 CO and [C I], we find support for the 
high-temperat ure molecular gas com ponent presented in 
the results of IPanuzzo et all (p OlO). The temperature 
traced by the warm component of 12 CO is quite high 
(335-518 K), and the addition of 13 CO slightly expands 
the likelihood ranges. The addition of [C I] produced 
results that indicate that these atom is not entirely trac- 
ing the same region as 12 CO. Some of the emission from 
these molecules (especially [C I]) are likely tracing more 
diffuse gas less shielded from UV radiation. 

The mapping observations did not resolve any signif- 
icant gradients in physical parameters (except evidence 
for a slight drop-off in beam-averaged column density, 
consistent with the beam profile) indicating that the sin- 
gle "deep" spectrum is an adequate representation of the 
galaxy when limited by our beam size. However, the 
mapping observations were important in confirming the 
source -beam coupling factor utilized in IPanuzzo et al.1 
(2010) and here, because through convolution of the 
maps we were able to confirm the central pixel's results 
matched with the deep spectrum. 

Molecular absorption traces lower column regions of 
the disk than those traced by CO emission, but con- 
tribute to the interpretation of the molecular gas of M82 
being excited by a combination of sources. Despite the 
enhanced cosmic ray density in M82, we do not find ev- 
idence that cosmic rays alone are sufficient to heat the 
gas enough to match the modeled hydrogen cooling rate. 
PDR models can only replicate the high-J CO line emis- 
sion at high densities incompatible with those indicated 
by the likelihood analysis, though cosmic-ray enhanced 
PDRs may be a closer match at lower densities. Turbu- 
lent heating from stellar winds and supernovae likely play 
a large role in the heating. More specifically, shocks are 
requi red to explain bright high-J line emission (|Pon et al.l 
l20ll . 

Like other submillimeter bright galaxies, Herschel has 
opened up new opportunities and questions about molec- 
ular and atomic lines that have never been observed be- 
fore. Because of this, the diagnostic power of high-J CO 
lines is still in development, and newer models currently 
being developed may be able to explain the emission seen 
from M82 and other extreme environments. 
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Fig. 2.1. — Integrated Flux (top) and Signal/Noise (bottom) maps for CI J = l— >0. 
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Fig. 2.2. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 5— >4. 
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Fig. 2.3. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 6— >5. 
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Fig. 2.4. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 7— s-6. 




Fig. 2.5. — Integrated Flux (top) and Signal/Noise (bottom) maps for CI J=2— 
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Fig. 2.6. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 8^7. 




Fig. 2.7. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 9— >8. 
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Fig. 2.8. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 10— >9. 
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Fig. 2.9. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 11— >10. 




Fig. 2.10. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 12— ►11. 
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Fig. 2.11. — Integrated Flux (top) and Signal/Noise (bottom) maps for NIL 
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Fig. 2.12. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 13— >12. 
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Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 4— >3. 
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Fig. 2.14. — Integrated Flux (top) and Signal/Noise (bottom) maps for CI J = l— s-0. 
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Fig. 2.15. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 5— v4. 



34 



KAMENETZKY ET AL. 



Arcseconds 




149.02 149.00 148.98 148.96 148.94 148.92 148.90 
RA 



Arcseconds 




Fig. 2.16. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 6— >5. 
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Fig. 2.17. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 7— >6. 
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Fig. 2.18. — Integrated Flux (top) and Signal/Noise (bottom) maps for CI J = 2— 
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Fig. 2.19. — Integrated Flux (top) and Signal/Noise (bottom) maps for CO J = 8— >7. 



